#! usr/bin/env python

"""We'll keep it simple for now. Takes a sequence and a codon usage table and optimizes it

Usage:

	./OptCodon.py input_file.txt

Arguments:
	input_file.txt ---> the file that specifies other files
"""
import xlwt
import random
import numpy as np
import re
import sys
import csv
import string
import math
import bisect
from string import maketrans
from sys import argv

__author__ = "Justin Finkle"
__credits__ = ["Justin Finkle"]
__version__ = "1.0"
__maintainer__ = "TBD"
__email__ = "jfinkle@u.northwestern.edu"
__status__ = "Prototype"


# Functions Functions Functions
#==============================================================================
#Split a string of characters into 3 base codons
def codon_splitter(sequence, chunk_size):
   return [sequence[i:i+chunk_size] for i in range(0, len(sequence), chunk_size)]

#Open the input file and get file locations, as a function for practice
def get_file_info(filename):
	try:
		infile = open(filename,'r')
	except IOError, e:
		sys.exit("No file by the name %s was found\nExiting..." %filename)
	file_info = infile.readlines()
	try:
		#Splits the input text file into parts
		if "Original_Sequence_File" not in file_info[0]:
			sys.exit("File location info not in correct format A")
		for entry in range(len(file_info)):
			divided_entry = file_info[entry].split("|")
			file_info[entry] = divided_entry[1].rstrip("\n")
		#Reads the information in each of these parts
		o_seq_file_location = file_info[0]
		codon_file_location=  file_info[1]
		gene_name = file_info[2]
		rare_codon_range = file_info[3]
		out_seq_file = file_info[4]
		excel_file_location = file_info[5]
	except Exception, e:
		sys.exit("File location info not in correct format B")
	return o_seq_file_location, codon_file_location, gene_name, rare_codon_range, out_seq_file, excel_file_location

def cdf_maker(probabilities):
	result = []  #Initialize blank list
	result = np.cumsum(probabilities)  #Populate with cumulative probabilities
	return result

#Pick a random codon based on probability distribution
def codon_selector(codons, probabilities):
    assert len(codons) == len(probabilities)  #Make sure they are the same length
    cdf_list=cdf_maker(probabilities)   #Get a CDF based on the PDF
    the_decision_maker=random.random()  #Random number
    
    #Using the random number, bisect the CDF and return the list index
    #Example: [0.1, 0.3, 0.7, 0.9, 1] with random number 0.8 returns 2
    index=bisect.bisect(cdf_list,the_decision_maker)
    if index>(len(cdf_list)-1):  #Sometimes the CDF doesn't add to 1. Slight to the last value
    	index=len(cdf_list)-1
    return codons[index]	
#End of Functions===============================================================

# Step 1a: Get the input file, get other file location and extract the original sequence
#==============================================================================
#Get the inputfile
script, input_filename = argv

#Get file info and open the file
o_seq_file, codon_file, gene_name, sweep_range,	output_seq_file, excel_file =get_file_info(input_filename)

#Get codon table name for use later
codon_table=codon_file.split('/')[-1].split('.')[0]

try:
	o_seq_infile = open(o_seq_file,'r')
except IOError, e:
	sys.exit("No sequence file by the name %s was found" %o_seq_file)

#Read the sequence into a list
o_seq_list = o_seq_infile.readlines()
o_seq_clean_list=[]
#Remove line breaks if they exist
for lines in o_seq_list:
	o_seq_clean_list.append(lines.strip())
#Join multiple lines
o_seq=''.join(o_seq_clean_list)  #Make it into a long string
#==============================================================================


# Step 1b: Basic Sequence debugging checks
#==============================================================================
#Verify that the sequence only contains the correct characters
allowed_chars = 'AaTtCcGg'
trans_table = string.maketrans('','')
invalid_chars_exist=o_seq.translate(trans_table, allowed_chars)

#If there are invalid characters, exit the script
if invalid_chars_exist:
	sys.exit('\nInvalid sequence characters\nExiting...')

#Get the length
o_seq_length = len(o_seq)
print ('The sequence entered has %d bases' %o_seq_length)

#Make sure it is the right length for codons
o_seq_codon_count = o_seq_length/3
extra_bases = o_seq_length%3

if (extra_bases>0):
	sys.exit("The sequence is not divisible by 3\nExiting...")

if (o_seq_codon_count>1):
	print ('The sequence entered has %d codons' %o_seq_codon_count)
else:
	print ('The sequence entered has 1 codon')
	
#Capitalize the sequence for ease of use
o_seq=string.upper(o_seq)
#==============================================================================


# Step 2: Read in a codon usage table file
#==============================================================================
codon_dict={}
master_dict={}
aa_dict={}
try:
	with open(codon_file, 'rU') as csvfile:  #variable codon_file is specified in step 1a
		codonreader = csv.DictReader(csvfile)
		for rows in codonreader:
			#format is: {Codon: {AA:,DNA:,Codon:,Number:,frequency:,fraction:}}
			master_dict[rows['DNA']] = rows  #Create a dictionary of dictionaries
			codon_dict[rows['DNA']] = 0 		        #Create a dictionary for codons to use later
			aa_dict[rows['AA']]=0						#Create a dictionary of aa for later
except IOError, e:
    sys.exit("No file by the name %s was found\nExiting..." %codon_file)

#Make a dictionary for the desired fraction
desired_codon_frac_dict={}
for codon, dict in master_dict.iteritems():
	desired_codon_frac_dict[codon] = float(dict['frac'])
	
#Make rare codon dictionary
o_rare_codon_count_dict={}
for codon, frac in desired_codon_frac_dict.iteritems():
	if frac==0:
		o_rare_codon_count_dict[codon]=0
#==============================================================================


#Step 3: Split the sequence into codons, get a count, and a frequency to compare
#		 Also, make a lot of different data structures, that may or may not be used
#==============================================================================
o_codon_list = codon_splitter(o_seq, 3)
o_codon_count_dict = codon_dict.copy()  #create a codon count dictionary

#Get original codon counts
for ii in o_codon_list:
	o_codon_count_dict[ii]+=1
	if ii in o_rare_codon_count_dict:
		o_rare_codon_count_dict[ii]+=1
		
#Calculate initial frequencies, per thousand		
o_codon_freq_dict = o_codon_count_dict.copy()
for codon,count in o_codon_freq_dict.iteritems():
	o_codon_freq_dict[codon]=round((o_codon_freq_dict[codon])/(o_seq_codon_count/1000.0),2)

#Make an amino acid library
codon_to_aa_dict={}
for codon, dict in master_dict.iteritems():
	if dict['DNA']==codon:
		codon_to_aa_dict[codon]=dict['AA']
#Sort it for future use

#Make an original AA list
o_aa_list=[]
for each_codon in range(len(o_codon_list)):
	for codon,aa in codon_to_aa_dict.iteritems():
		if o_codon_list[each_codon]==codon:
			o_aa_list.append(aa)

#Get a count of the original AA
o_aa_count_dict = aa_dict.copy()  #initialize it to a blank aa dict
for aa, count in o_aa_count_dict.iteritems():
	for o_aa in o_aa_list:
		if o_aa==aa:
			o_aa_count_dict[aa]+=1

#Calculate original fraction for reporting
o_codon_frac_dict = o_codon_count_dict.copy()		#Make an empty codon dict
for codon, fraction in o_codon_frac_dict.iteritems():
	#Catch if count is 0 fraction division won't work
	if o_aa_count_dict[codon_to_aa_dict[codon]]==0:
		o_codon_frac_dict[codon]=0
	else:
		o_codon_frac_dict[codon] = round(o_codon_frac_dict[codon]/float(o_aa_count_dict[codon_to_aa_dict[codon]]),2)
#print o_codon_frac_dict
#==============================================================================


#Step 4: Get stats to decide if optimization is worth it
#==============================================================================
o_codon_count = sum(o_codon_count_dict.values())
o_rare_codon_count = sum(o_rare_codon_count_dict.values())
if o_codon_count!=o_seq_codon_count:
	sys.exit("Number of codons does not match expected\nExiting...")
o_rare_percent = round(o_rare_codon_count/float(o_codon_count),4)*100

#Get the original positions of the rare codons
o_rare_positions=[]
for codon_index in range(len(o_codon_list)):
	if o_codon_list[codon_index] in o_rare_codon_count_dict:
		o_rare_positions.append(codon_index)

close_rare_count=len(o_rare_positions)
#Make sure sweep is an integer
sweep_range=int(sweep_range)
close_rare_dict = {}
temp_list =[]
for ii in range(len(o_rare_positions)-1):
	base_position = o_rare_positions[ii]
	next_position = o_rare_positions[ii+1]
	if (next_position-base_position)<=sweep_range:
		temp_list.append(base_position)
	else:
		if (base_position-o_rare_positions[ii-1])<=sweep_range:
			temp_list.append(base_position)
		if len(temp_list)>=1:
			close_rare_dict[temp_list[0]]=temp_list[:]
		temp_list=[]

#Get lengths of stretches for later
rare_print_info={}
for start, list in close_rare_dict.iteritems():
	list_count=len(list)
	if list_count in rare_print_info:
		rare_print_info[list_count]+=1
	else:
		rare_print_info[list_count]=1
sorted_rare_print_list=sorted(rare_print_info.iterkeys())
#==============================================================================

#Step 5: Create a dictionary of dictionaries that contains the data we will use
#==============================================================================

#Create a dictionary that has aa pertinent info
o_aa_info_dict=aa_dict.copy()	#initialize a blank aa dict to hold our info

#NOTE: dict.fromkeys() did not seem to work. It appeared to link dictionaries from different
#codons together. 
for aa, info in o_aa_info_dict.iteritems():
	o_aa_info_dict[aa]={'associated_codons':[], 'desired_fraction':[],
						'count':0, 'ideal_count':[]}  #new codon will be populated later

#Make a comprehensive library for later use
for aa, info in o_aa_info_dict.iteritems():
	
	#Get a count to be safe
	for codon_pos in range(len(o_aa_list)):
		if o_aa_list[codon_pos]==aa:
			o_aa_info_dict[aa]['count']+=1
	#Get the right codons and fractions
	for lookup_codon, lookup_aa in codon_to_aa_dict.iteritems():
		if lookup_aa==aa:
			o_aa_info_dict[aa]['associated_codons'].append(lookup_codon)
			#Catch if count is 0 fraction division won't work
			if o_aa_info_dict[aa]['count']==0:
				rounding_factor=0
			else:
				rounding_factor = float(1)/o_aa_info_dict[aa]['count']
			for desired_codon, fraction in desired_codon_frac_dict.iteritems():
				if desired_codon==lookup_codon:
					float_fraction=float(fraction)
					o_aa_info_dict[aa]['desired_fraction'].append(float_fraction) #need to make them floats
					
					#Round to the nearest (1/sequence_aa_count). This is the the best
					#fraction that can be obtained for the specific sequence given
					#Catch if count is 0 fraction division won't work
					if rounding_factor==0:
						ideal_count=0
					else:
						ideal_count = int(round(float_fraction/rounding_factor,0))
					o_aa_info_dict[aa]['ideal_count'].append(ideal_count)
	
	#Distribute any rounding errors accordingly
	ideal_sum = sum(o_aa_info_dict[aa]['ideal_count'])
	real_sum = o_aa_info_dict[aa]['count']
	while ideal_sum!=real_sum:
		pick_list=[]
		for index in range(len(o_aa_info_dict[aa]['desired_fraction'])):
			if o_aa_info_dict[aa]['desired_fraction'][index]>0:
				pick_list.append(index)
		random_select = random.choice(pick_list)
		if ideal_sum>real_sum:
			o_aa_info_dict[aa]['ideal_count'][random_select]+=(-1)
		else:
			o_aa_info_dict[aa]['ideal_count'][random_select]+=1
		#Recalc sums
		ideal_sum = sum(o_aa_info_dict[aa]['ideal_count'])
		real_sum = o_aa_info_dict[aa]['count']
		
#Build the ideal count dictionary
ideal_count_dict={}
for aa, info in o_aa_info_dict.iteritems():
	for codon_index in range(len(o_aa_info_dict[aa]['associated_codons'])):
		ideal_count_dict[o_aa_info_dict[aa]['associated_codons'][codon_index]] \
		=o_aa_info_dict[aa]['ideal_count'][codon_index]
#==============================================================================


#Step 6: Optimize. At least, that's the plan...
#==============================================================================

#Create copies to be safe
opt_aa_info_dict=o_aa_info_dict.copy()
opt_aa_list=o_aa_list[:]
opt_codon_list=o_codon_list[:]
opt_codon_frac_dict = dict.fromkeys(o_codon_count_dict.keys(),0)
opt_codon_count_dict = o_codon_count_dict.copy()
opt_difference_dict = dict.fromkeys(opt_codon_count_dict.keys(),0)

#Make a blank list of indices for reference
index_list=[]
for ii in range(0,len(opt_aa_list)):
	index_list.append(ii)

loop_count=0
changes=1
totaldiff=1
old_totaldiff=0
#while (old_totaldiff!=totaldiff) and (totaldiff>0.06):
#for jj in range(1,1000):
while (opt_codon_count_dict!=ideal_count_dict):
	#Randomize index list for selection
	loop_count+=1
	print 'Loop #: %d' %loop_count
	rand_index_list=random.sample(index_list, len(index_list))
	changes=0
	for kk in range(len(rand_index_list)):
		#Get the current sequence index from the randomized list
		current_seq_index=rand_index_list[kk]

		#Get the current aa, codon, and fraction
		current_aa=opt_aa_list[current_seq_index]
		current_codon= opt_codon_list[current_seq_index]
		current_count = opt_codon_count_dict[current_codon]
		
		#Get the current aa dictionary for ease of use
		current_aa_dict=opt_aa_info_dict[current_aa].copy()
		
		#Get the desired optimal fraction
		opt_count_index = current_aa_dict['associated_codons'].index(current_codon)
		opt_count = current_aa_dict['ideal_count'][opt_count_index]
		
# 		print 'current_index: %d' %current_seq_index
# 		print 'current_aa: '+ current_aa
# 		print 'current_codon: ' + current_codon
# 		print 'current_count: %f' %current_count
# 		print current_aa_dict
# 		print opt_count
# 		raw_input()
		
		#If the current count is greater than the desired count, optimize!
		if current_count > opt_count:
			#print 'Optimize!'
			#Get the desired counts to use as a probability
			codon_counts=current_aa_dict['ideal_count'][:]
			
			#Normalize the counts for probability selection
			#This could be done once before hand to improve speed
			n_fractions=[] #Reinitialize
			for each_item in range(len(codon_counts)):
				n_fractions.append(codon_counts[each_item]/float(sum(codon_counts)))
# 			print n_fractions
# 			print sum(n_fractions)
# 			raw_input()
			available_codons=current_aa_dict['associated_codons']
			#print fractions
			#print available_codons
			new_codon=codon_selector(available_codons,n_fractions)
			#print new_codon
			opt_codon_list[current_seq_index]=new_codon
			if new_codon!=current_codon:
				#Log the change
				changes+=1
		
		#Recalculate fractions
		
		#Recount each codon
		opt_codon_count_dict = dict.fromkeys(ideal_count_dict.keys(),0) #Reinitialize
		opt_aa_count_dict=dict.fromkeys(o_aa_count_dict.keys(),0)
		for codon,count in opt_codon_count_dict.iteritems():
			for i in opt_codon_list:
				if i==codon:
					opt_codon_count_dict[codon]+=1
					opt_aa_count_dict[codon_to_aa_dict[codon]]+=1
	
	print 'Changes: %d'%changes
	
	#Calculate count difference
	old_totaldiff=totaldiff
	totaldiff=0
	for opt_codon,opt_count in opt_codon_count_dict.iteritems():
		#print opt_codon]
		#print opt_fraction
		#print desired_opt_frac_dict[opt_codon]
		#raw_input()
		opt_difference_dict[opt_codon]=opt_count-ideal_count_dict[opt_codon]
		totaldiff=abs(opt_difference_dict[opt_codon])+totaldiff
		#print opt_difference_dict[opt_codon]
		#raw_input()
	#totaldiff=sum(opt_difference_dict.itervalues())
	print 'Differences: %d' %totaldiff
	#print opt_codon_count_dict
	#print ideal_count_dict
	#print opt_difference_dict
	print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
#Calculate the new fraction
for codon, fraction in opt_codon_frac_dict.iteritems():
	#Catch if count is 0 fraction division won't work
	if opt_aa_count_dict[codon_to_aa_dict[codon]]==0:
		opt_codon_frac_dict[codon]=0
	else:
		opt_codon_frac_dict[codon] = round(opt_codon_count_dict[codon]/float(opt_aa_count_dict[codon_to_aa_dict[codon]]),2)

#print 'Differences: %d' %totaldiff
#==============================================================================


#Step 7: Cleanup and checks
#==============================================================================
print_to_file=0 #This changes to one if the sequence passes the subsequent tests
print 'Optimization run details:'
print '\tIterations Taken: %d' %loop_count

#Make an optimized aa list
opt_aa_list_check=[]
for each_codon in range(len(opt_codon_list)):
	for codon,aa in codon_to_aa_dict.iteritems():
		if opt_codon_list[each_codon]==codon:
			opt_aa_list_check.append(aa)

#Get a count of the optimized AA
opt_aa_count_dict_check = dict.fromkeys(opt_aa_count_dict.keys(),0)  #initialize it to a blank aa dict
for aa, count in opt_aa_count_dict_check.iteritems():
	for opt_aa in opt_aa_list_check:
		if opt_aa==aa:
			opt_aa_count_dict_check[aa]+=1
			
#Verify the count
if opt_aa_count_dict_check==o_aa_count_dict:
	print '\tAA count dictionaries match'
	#Verify the list
	aa_diff=0
	aa_diff_loc=[]
	for kk, aa in enumerate(opt_aa_list_check):
		if aa==o_aa_list[kk]:
			pass
		else:
			aa_diff+=1
			aa_diff_loc.append(kk)
	if aa_diff==0:
		print '\tAA sequences match'
		print_to_file=1
		#Get a string list
		opt_seq_string=''.join(opt_codon_list)
	else:
		print '\tAA sequence differ at these locations %r' %aa_diff_loc
else:
	print '\tAA count dictionaries do not match'

same_bases=0
for ii,base in enumerate(o_seq):
	if base==opt_seq_string[ii]:
		pass #print base+'==='+opt_seq_string[i]
	else:
		#print base+'XXX'+opt_seq_string[i]
		same_bases+=1
print '\tUnchanged Bases: %d' %same_bases
#==============================================================================


#Step 7: Print to file
#==============================================================================
if print_to_file==1:
	print 'Writing to file...'
	
	#Write the optimized sequence to a txt file
	seq_write_file = open(output_seq_file, 'w')
	seq_write_file.write(opt_seq_string)
	seq_write_file.close()

	#Initialize workbook
	wbk = xlwt.Workbook()
	
	#Make necessary sheet(s)
	results_sheet = wbk.add_sheet('Optimizer Results')
	codon_sheet = wbk.add_sheet('Codon Comparison')
	sequence_sheet = wbk.add_sheet('Sequence Comparison')
	
	#Create styles
	style_default=xlwt.easyxf('alignment: horizontal center;')
	style_title=xlwt.easyxf('font: bold 1;' 'alignment: horizontal center, vertical center, wrap 1;')
	style_green=xlwt.easyxf('font: color green;' 'alignment: horizontal center;')
	style_red=xlwt.easyxf('font: color red;' 'alignment: horizontal center;')
	style_bold=xlwt.easyxf('font: bold 1;')
	style_rare =xlwt.easyxf('pattern: pattern solid, fore_color red;' 'alignment: horizontal center;')
	style_bold_left=xlwt.easyxf('font: bold 1;' 'alignment: horizontal left;')
	style_left=xlwt.easyxf('alignment: horizontal left;')
	#Height of 320 = 16 pt font, 1pt=20 for some reason
	style_big_title=xlwt.easyxf('font: bold 1, height 320;' 'alignment: horizontal left;'
								'border: bottom double;')
#=======================Write sequence info to file==============================	
	#Set current sheet
	current_sheet = sequence_sheet
	#Write sequence info to current_sheet
	#>>>>>>>>>>>>>Write format is sheet.write(row, column, value, style)<<<<<<<<<<<<<<<<<
	#Write titles
	current_sheet.write(0,0,'Base #', style_title)
	current_sheet.write(0,1,'Original Sequence', style_title)
	current_sheet.write(0,3,'Optimized Sequence', style_title)

	#NOTE: The first row is used for headers, will need to offset by 1
	for ii,base in enumerate(o_seq):
		current_sheet.write(ii+1,0, ii+1, style_default) 					#Base number
		current_sheet.write(ii+1,1, o_seq[ii], style_default) 				#Original base
		current_sheet.write(ii+1,3, opt_seq_string[ii], style_default) 		#Optimized base
		#Print whether the original and optimized base are the same or not
		if base==opt_seq_string[ii]:
			current_sheet.write(ii+1,2, '===', style_green)
		else:
			current_sheet.write(ii+1,2, '=/=', style_red)
#==============================================================================	

#=======================Write codon sequence info to file==============================
	#Set current sheet
	current_sheet = codon_sheet
	
	#Write titles
	current_sheet.write(0,0,'Codon #', style_title)
	current_sheet.write(0,1,'Original Codons', style_title)
	current_sheet.write(0,3,'Optimized Codons', style_title)
	
	codons_changed=0  #Initialize codon change counter
	for ii,codon in enumerate(o_codon_list):
		#If it's a rare codon, give it a red background
		current_sheet.write(ii+1,0, ii+1, style_default)				#Base number
		#Print highlighted cell if it is originally a rare codon
		if ii in o_rare_positions:
			current_sheet.write(ii+1,1, o_codon_list[ii], style_rare)
		else:
			current_sheet.write(ii+1,1, o_codon_list[ii], style_default)
		#Print the new codon	
		current_sheet.write(ii+1,3, opt_codon_list[ii], style_default)
		if codon==opt_codon_list[ii]:
			current_sheet.write(ii+1,2, '===', style_green)
		else:
			current_sheet.write(ii+1,2, '=/=', style_red)
			codons_changed+=1   #Tally a new codon
#===================================================================================

#========================Write codon stat info to file==============================	
	#Set current sheet
	current_sheet =results_sheet
		
	#Statistics offset for the number of rows
	stat_offset = 6
	
	titles_list = ['Codon','AA', 'Original Count', 'Original Fraction', \
					'Desired Fraction', 'Optimized Count', 'Optimized Fraction']
					
	for ii in range(len(titles_list)):
		current_sheet.write(stat_offset, ii, titles_list[ii], style_title)	
	
	write_row=0
	sorted_aa_list=sorted(codon_to_aa_dict.iterkeys())
	for codon in sorted_aa_list:
		current_sheet.write(write_row+stat_offset+1, 0, codon, style_default)
		current_sheet.write(write_row+stat_offset+1, 1, codon_to_aa_dict[codon], style_default)
		current_sheet.write(write_row+stat_offset+1, 2, o_codon_count_dict[codon], style_default)
		current_sheet.write(write_row+stat_offset+1, 3, o_codon_frac_dict[codon], style_default)
		current_sheet.write(write_row+stat_offset+1, 4, desired_codon_frac_dict[codon], style_default)
		current_sheet.write(write_row+stat_offset+1, 5, opt_codon_count_dict[codon], style_default)
		current_sheet.write(write_row+stat_offset+1, 6, opt_codon_frac_dict[codon], style_default)
		write_row+=1
	
	#Write Statistics
	#Header
	header_list = ['Optimizer Statistics', '', '', '', '','']
	for ii in range(len(header_list)):
		current_sheet.write(0, ii, header_list[ii], style_big_title)
		
	#Write misc info
	percent_changed = codons_changed*100/o_seq_codon_count
	percent_rare = close_rare_count*100/o_seq_codon_count
	
	current_sheet.write(1,0, 'Gene Name', style_bold_left)
	current_sheet.write(1,1, gene_name, style_left)		
	current_sheet.write(2,0, 'Seq Bases', style_bold_left)
	current_sheet.write(2,1, o_seq_length, style_left)
	current_sheet.write(3,0, 'Seq Codons', style_bold_left)
	current_sheet.write(3,1, o_seq_codon_count, style_left)
	
	current_sheet.write(1,3, 'Optimization Loops', style_bold_left)
	current_sheet.write(1,5, loop_count, style_left)
	current_sheet.write(2,3, 'Codons Changed', style_bold_left)
	current_sheet.write(2,5, '%r (%r%%)' %(codons_changed, percent_changed), style_left)
	current_sheet.write(3,3, 'Rare codons', style_bold_left)
	current_sheet.write(3,5, '%r (%r%%)' %(close_rare_count, percent_rare), style_left)
	current_sheet.write(4,3, 'Optimized using', style_bold_left)
	current_sheet.write(4,5, codon_table, style_left)
	
	#Write rare info
	if close_rare_count>0:
		gap = len(titles_list)+1
		rare_header_list = ['Rare Cluster Length', '# in Sequence']
		for ii in range(len(rare_header_list)):
			current_sheet.write(stat_offset, gap+ii, rare_header_list[ii], style_title)
		for kk in range(len(sorted_rare_print_list)):
			current_sheet.write(stat_offset+1+kk,gap+0, sorted_rare_print_list[kk], style_default)
			current_sheet.write(stat_offset+1+kk,gap+1, rare_print_info[sorted_rare_print_list[kk]], style_default)
	
	wbk.save(excel_file)
	print 'File wrote to location %r' %excel_file
	print 'Optimization Complete!'
